Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation.

This script was used to undertake the analysis of the primary outcome for ASCOT.

Author

James Totterdell

Published

July 14, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)
library(ggdist)
library(patchwork)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))
color_scheme_set("red")
options(digits = 4)
Prepare datasets
devtools::load_all()
# Raw data
all_dat <- read_all_no_daily()
all_dat_daily <- read_all_daily()

# FAS-ITT
fas_itt_dat <- all_dat %>% 
  filter_fas_itt() %>%
  transmute_model_cols_grp_aus_nz()
# FAS-ITT excluding missing values
fas_itt_nona_dat <- fas_itt_dat %>% 
  filter(!is.na(PO))

# ACS-ITT
acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()
# ACS-ITT excluding missing values
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(PO))
# ACS-ITT worst case
acs_itt_wc_dat <- acs_itt_dat %>%
  mutate(PO = if_else(is.na(PO) | PO == 1, 1, 0))
# ACS-ITT best case
acs_itt_bc_dat <- acs_itt_dat %>%
  mutate(PO = if_else(is.na(PO) | PO == 0, 0, 1))

# ACS-PP - exclude participants not per-protocol
pp <- read_per_protocol_list()
all_pp <- all_dat %>% left_join(pp, by = "StudyPatientID")
acs_pp_dat <- all_pp %>%
  filter_acs_itt() %>%
  filter(is.na(Reason)) %>%
  transmute_model_cols_grp_aus_nz()
acs_pp_nona_dat <- acs_pp_dat %>%
  filter(!is.na(PO))

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(PO))
# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(PO))
# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(PO))
Load models
logistic_mod <- cmdstan_model("stan/binary/logistic.stan")
logistic_site_mod <- cmdstan_model("stan/binary/logistic_site.stan")
logistic_epoch_mod <- cmdstan_model("stan/binary/logistic_epoch.stan")
logistic_site_epoch_llik <- cmdstan_model("stan/binary/logistic_site_epoch_loglik.stan")
logistic_site_epoch_llik_int <- cmdstan_model("stan/subgroup/logistic_site_epoch_loglik_shrinkage.stan")
primary_mod <- cmdstan_model("stan/binary/logistic_site_epoch.stan")
Helper functions
make_stan_data <- function(
    dat, 
    vars, 
    outcome, 
    beta_sd = c(2.5, rep(1, nXassign), 2.5, 10),
    ctr = contr.orthonorm) 
{
  X <- model.matrix(
     as.formula(paste("~", paste(vars, collapse = " + "))),
    data = dat,
    contrast = list(
      AAssignment = ctr,
      CAssignment = ctr)
  )
  nXassign <- sum(grepl("Assignment", colnames(X)))

  y <- dat[[outcome]]
  N <- dim(X)[1]
  K <- dim(X)[2]  
    epoch  <- dat$epoch
  M_epoch  <- max(dat$epoch)
    region <- dat[["ctry_num"]]
  M_region <- max(region)
    site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  list(N = N, K = K, X = X, y = y,
       M_region = M_region, region = region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       beta_sd = beta_sd)
}


# Always includes the assigned interventions
make_primary_model_data <- function(
    dat, 
    vars = c("inelgc3", "agegte60", "ctry"),
    beta_sd_var = c(10, 2.5, 1, 1),
    ctr = contr.orthonorm) {
  
  # Domain specific intercept as in interims
  XA <- make_domA_design(dat, ctr)
  # Domain specific intercept as in interims
  XC <- make_domC_design(dat, ctr)
  if (any(XC[, 1] != 1)) {
    add_intercept = TRUE
  } else {
    add_intercept = FALSE
  }
  
  if (!is.null(vars)) {
    Xother <- model.matrix(
     as.formula(paste(" ~ ", paste(vars, collapse = " + "))),
     data = dat
    ) [, -1]
  } else {
    Xother <- NULL
  }
  
  X <- cbind(XC, XA, Xother)
  if(add_intercept) {
    X <- cbind(randPlatform = 1, X)
  }
  nXtrt <- sum(grepl("rand", colnames(X))) - 1
  
  epoch  <- dat[["epoch"]]
  M_epoch  <- max(epoch)
  region <- dat[["ctry_num"]]
  M_region <- max(region)
  site <- dat[["site_num"]]
  M_site <- max(site)
  region_by_site <- region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)
  
  y <- dat[["PO"]]
  N <- dim(X)[1]
  K <- dim(X)[2]
  beta_sd <- c(2.5, rep(1, nXtrt), beta_sd_var)
  out <- list(
    N = N, K = K, X = X, y = y,
    M_region = M_region, region = region,
    M_site = M_site, site = site,
    M_epoch = M_epoch, epoch = epoch,
    region_by_site = region_by_site,
    beta_sd = beta_sd,
    ctrC = attr(XC, "contrasts")[[1]],
    ctrA = attr(XA, "contrasts")[[1]]
  )
  return(out)
}


fit_primary_model <- function(
    dat,
    model = primary_mod,
    vars = NULL,
    beta_sd_var = NULL,
    ctr = contr.orthonorm,
    seed = 32915,
    ...
) {
  mdat <- make_primary_model_data(
    dat,
    vars,
    beta_sd_var,
    ctr
  )
  snk <- capture.output(
    mfit <- model[["sample"]](
      data = mdat,
      seed = seed,
      refresh = 0,
      iter_warmup = 1000,
      iter_sampling = 2500,
      chains = 8,
      parallel_chains = min(8, parallel::detectCores()),
      ...
  ))
  mpars <- mfit$metadata()$model_params
  keep <- mpars[!grepl("(_raw|epsilon_)", mpars)]
  mdrws <- as_draws_rvars(mfit$draws(keep))
  names(mdrws$beta) <- colnames(mdat$X)
  
  if(any(grepl("site", keep))) {
    site_map <- dat %>% dplyr::count(site_num, site)  
    names(mdrws$gamma_site) <- site_map$site
  }
  if(any(grepl("epoch", keep))) {
    epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)  
    names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  }
  # Transformed samples
  Ca <- mdat$ctrA
  Cc <- mdat$ctrC
  # Get constrained parameters from contrast transformation
  mdrws$Acon <- Ca %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- as.vector(Cc %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))])
  names(mdrws$Ccon) <- rownames(Cc)
  # Get treatment effect in case of non treatment-coding
  mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  # Get odds ratios
  mdrws$OR <- exp(mdrws$Ctrt)
  mdrws$OR <- c(mdrws$OR, exp(mdrws$beta[!grepl("(Intercept|rand)", names(mdrws$beta))]))
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}


make_po_table <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"), 
      labels = intervention_labels()$CAssignment[-1])
    ) %>%
    group_by(CAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(PO)), 100 * mean(is.na(PO))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
      `Met primary outcome` = sprintf("%i (%.1f)", sum(PO, na.rm = TRUE), 100 * mean(PO, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(PO)), 100 * mean(is.na(PO))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
          `Met primary outcome` = sprintf("%i (%.1f)", sum(PO, na.rm = TRUE), 100 * mean(PO, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(CAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}


summarise_posterior <- function(dat, futval = 1.1) {
  dat %>%
    mutate(
      Median = 
        sprintf("%.2f", median(Posterior)),
      `95\\% CrI` = 
        sprintf("(%.2f, %.2f)", quantile(Posterior, 0.025), quantile(Posterior, 0.975)),
      `Mean (SD)` = 
        sprintf("%.2f (%.2f)", mean(Posterior), sd(Posterior)),
      `Pr(OR < 1)` = 
        sprintf("%.2f", Pr(Posterior < 1)),
      `Pr(OR > 1/1.1)` = 
        sprintf("%.2f", Pr(Posterior > 1/futval))
    )
}


make_subgroup_summary_table <- function(dat, format = "html") {
  L <- nrow(dat)
  row_st <- seq(1, L, by = L / 3)
  row_ed <- seq(L/3, L, by = L / 3)
  kable(
    dat[, -1] %>% select(-Posterior),
    format = format,
    booktabs = TRUE,
    align = "lrrrrr",
      escape = FALSE
  ) %>%
    kable_styling(
      latex_options = "HOLD_position",
      font_size = 9
    ) %>%
    group_rows("Intermediate", row_st[1], row_ed[1]) %>%
    pack_rows("Low with aspirin", row_st[2], row_ed[2]) %>%
    pack_rows("Therapeutic", row_st[3], row_ed[3])  
}

plot_or_densities <- function(rvs) {
  tibble(Contrast = fct_inorder(names(rvs)), RV = rvs) %>% 
  ggplot(., aes(y = Contrast, xdist = RV)) + 
  stat_halfeye(
    aes(fill = 
          stat(cut_cdf_qi(
            cdf, 
            .width = c(.5, .8, .95, 0.99), 
            labels = scales::percent_format()))),
    adjust = 1, n = 1001, .width = c(0.5, 0.8, 0.95)
    ) + 
  scale_fill_brewer(
    palette = "Reds",
    direction = -1, 
    na.translate = FALSE) + 
  labs(
    x = "Odds ratio contrast",
    fill = "Interval"
  ) +
  scale_x_log10(breaks = c(0.1, 0.25, 0.5, 1, 2, 4, 10)) + 
  geom_vline(xintercept = 1)
}


odds_ratio_summary_table <- function(OR, format = "html", fn = NULL) {
  out <- tibble(
    Parameter = names(OR),
    Median = median(OR),
    `95% CrI` = apply(
      quantile(OR, c(0.025, 0.975)), 2, 
      function(x) sprintf("(%.2f, %.2f)", x[1], x[2])),
    `Mean (SD)` = sprintf("%.2f (%.2f)", E(OR), sd(OR)),
    `Pr(OR < 1)` = Pr(OR < 1),
  ) %>%
  kable(
    format = format, 
    digits = 2, 
    align = "lrrrr",
    linesep = "",
    booktabs = TRUE) %>%
  kable_styling(
    font_size = 9,
    bootstrap_options = "striped", 
    latex_options = "HOLD_position")
  if (!is.null(fn) & format == "latex") {
    save_tex_table(out, fn)
  } else {
    return(out)
  }
}


decision_quantity_summary_table <- function(OR, format = "html", fut_val = 1.1) {
  tdat <- tibble(
    Intervention = names(OR),
    posterior = OR,
    Posterior = sprintf("%.2f (%.2f, %.2f)", 
                        median(posterior), 
                        quantile(posterior, 0.025), 
                        quantile(posterior, 0.975)),
    `Superior\nPr(OR = min(OR))` = sprintf("%.2f", E(posterior == rvar_min(posterior))),
    `Effective\nPr(OR < 1)` = sprintf("%.2f", Pr(posterior < 1)),
    `Futile\nPr(OR > 1/1.1)` = sprintf("%.2f", Pr(posterior > 1/fut_val)),
    `Equivalent\nPr(1/1.1 < OR < 1.1)` = sprintf("%.2f", Pr(posterior < fut_val & posterior > 1/fut_val))
  ) %>%
    select(-posterior)
  tdat[1, 2] <- "1.00"
  tdat[1, -(1:3)] <- "-"
  if(format == "latex") {
    colnames(tdat) <- linebreak(colnames(tdat), align = "c", linebreaker = "\n")  
  }
  out <- kable(
    tdat,
    format = format, 
    digits = 2, 
    align = "lrrrrr",
    linesep = "",
    escape = F,
    booktabs = TRUE) %>%
    kable_styling(
      font_size = 9,
      bootstrap_options = "striped",
      latex_options = "HOLD_position"
    )
  return(out)
}

Outline

This document presents the results of analyses of the primary outcome. The first section introduces the primary outcome definition and it’s derivation from the database fields. Descriptive analyses are then presented followed by the model-based analyses.

For the model-based analyses, results are presented for each population and sub-model.

Primary Outcome Definition

The primary outcome is a composite of death, need for new respiratory support, or vasopressor/inotropic support. From the protocol:

Death from any cause or requirement of new intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support in the 28 days after randomisation. This includes any participant who receives non-invasive mechanical ventilation (either CPAP or BIPAP, apart from the below considerations) any time after enrolment even if not transferred to ICU. It does NOT include the use of humidified high-flow nasal prong oxygen.

Participants on pre-existing home BiPAP or CPAP will not be considered to have met the primary outcome unless they have either:

  • required invasive mechanical ventilation (i.e. intubation), or
  • graduated from CPAP only whilst asleep to BiPAP at any time, or
  • graduated from BiPAP only whilst asleep to BiPAP for >12 hours/day, or
  • died by day 28

There may be cases where a patient has been assessed as requiring intensive respiratory support (invasive or non-invasive ventilation) or vasopressor/inotropic support, but the patient or family declined treatment and the patient was discharged home. If attempts to obtain 28-day data are unsuccessful or not possible, and the investigator had deemed at the time of discharge that the patient would be highly likely to die within 28 days from randomisation, these participants will be deemed to have met the primary outcome.

Derivation of the Outcome

Derivation of the outcome requires checking of the daily, discharge, and day 28 data extracts. On the daily data, there is a variable DD_PrimaryEndpointReachedToday however this was coded incorrectly in the original database and therefore fails to capture some participants. Additionally, given the composite nature of the outcome it is useful to check all components individually as well as the composite outcome. Therefore, this variable is not used in the derivation, but is included as a cross-check.

Below, each component is summarised in aggregate.

Day 28 mortality

Summarise the death component
mort_table <- all_dat %>%
  filter_acs_itt() %>%
  group_by(CAssignment) %>%
  summarise(
    `Mortality_Alive at day 28_` = sum(1 - D28_death, na.rm = TRUE),
    `Mortality_Death within 28 days_Total` = sum(D28_death, na.rm = TRUE),
    `Mortality_Death within 28 days_Prior to discharge` = sum(DIS_death & DIS_day <= 28, na.rm = TRUE),
    `Mortality_Death within 28 days_Post-discharge` = sum(D28_death == 1 & DIS_death == 0, na.rm = TRUE),
    `Mortality_Unknown_` = sum(is.na(D28_death))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[c(1,2,5)])))
  )
kable(
  mort_table %>% select(-1),
  caption = "Composite primary endpoint breakdown - mortality.",
  align = "llrrrrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Mortality", 1, nrow(mort_table))
Outcome Breakdown C1 C2 C3 C4 Overall
Mortality
Alive at day 28 577 (95) 588 (96) 271 (96) 44 (88) 1480 (95)
Death within 28 days Total 19 (3) 15 (2) 10 (4) 6 (12) 50 (3)
Death within 28 days Prior to discharge 15 (2) 11 (2) 10 (4) 4 (8) 40 (3)
Death within 28 days Post-discharge 4 (1) 4 (1) 0 (0) 2 (4) 10 (1)
Unknown 14 (2) 10 (2) 2 (1) 0 (0) 26 (2)
Composite primary endpoint breakdown - mortality.

Vasopressor/inotropic requirement

Summarise the vasopressor/inotropes component
vaso_table <- all_dat %>%
  filter_acs_itt() %>%
  group_by(CAssignment) %>%
  summarise(
    `Vasopressor/inotropes_Not required_` = sum(1 - ANY_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Use within 28 days_Total` = sum(ANY_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Use within 28 days_Prior to discharge` = sum(DD_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Use within 28 days_Post-discharge` = sum(D28_vasop, na.rm = TRUE),
    `Vasopressor/inotropes_Unknown_` = sum(is.na(ANY_vasop))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[c(1,2,5)])))
  )
kable(
  vaso_table %>% select(-1),
  caption = "Composite primary endpoint breakdown - vasopressor/inotropes.",
  align = "llrrrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Vasopressor/inotropes", 1, nrow(vaso_table))
Outcome Breakdown C1 C2 C3 C4 Overall
Vasopressor/inotropes
Not required 591 (97) 595 (97) 273 (96) 48 (96) 1507 (97)
Use within 28 days Total 6 (1) 6 (1) 6 (2) 2 (4) 20 (1)
Use within 28 days Prior to discharge 5 (1) 6 (1) 6 (2) 2 (4) 19 (1)
Use within 28 days Post-discharge 2 (0) 0 (0) 0 (0) 1 (2) 3 (0)
Unknown 13 (2) 12 (2) 4 (1) 0 (0) 29 (2)
Composite primary endpoint breakdown - vasopressor/inotropes.

New intensive respiratory support

Summarise the ventilation component
vent_table <- all_dat %>%
  filter_acs_itt()  %>%
  group_by(CAssignment) %>%
  summarise(
    `Ventilation_Not required_` = sum(1 - ANY_vent, na.rm = TRUE),
    `Ventilation_Use within 28 days_Total` = sum(ANY_vent, na.rm = TRUE),
    `Ventilation_Use within 28 days_Prior to discharge` = sum(ANY_DD_vent, na.rm = TRUE),
    `Ventilation_Use within 28 days_Post-discharge` = sum(ANY_vent & ANY_DD_vent == 0, na.rm = TRUE),
    `Ventilation_Unknown_` = sum(is.na(ANY_vent))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[c(1, 2, 5)])))
  )
kable(
  vent_table %>% select(-1),
  caption = "Composite primary endpoint breakdown - respiratory support.",
  align = "llrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Ventilation", 1, nrow(vent_table))
Outcome Breakdown C1 C2 C3 C4 Overall
Ventilation
Not required 562 (92) 579 (94) 263 (93) 43 (86) 1447 (93)
Use within 28 days Total 35 (6) 24 (4) 18 (6) 7 (14) 84 (5)
Use within 28 days Prior to discharge 29 (5) 19 (3) 13 (5) 5 (10) 66 (4)
Use within 28 days Post-discharge 6 (1) 5 (1) 5 (2) 2 (4) 18 (1)
Unknown 13 (2) 10 (2) 2 (1) 0 (0) 25 (2)
Composite primary endpoint breakdown - respiratory support.

DAMA and Likely to Die

There were 3 participants who were DAMA and identified as likely to die. However, the day 28 status was observed for all 3 participants, therefore, no-one met the primary endpoint due to DAMA and unknown day 28 status.

Summarise the DAMA component
all_dat %>%
  filter_acs_itt()  %>%
  select(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
  dplyr::count(DIS_DAMA, DIS_DAMAlikelytodie, D28_death) %>%
  spread(D28_death, n) %>%
  kable(
    col.names = c("DAMA", "Likely to die by day 28", "No", "Yes", "(Missing)")
  ) %>%
  kable_styling() %>%
  add_header_above(c(" " = 2, "Death by day 28" = 3))
Death by day 28
DAMA Likely to die by day 28 No Yes (Missing)
0 0 1447 46 24
1 0 32 2 2
1 1 1 2 NA

Composite

Summarise the composite
comp_table <- all_dat %>%
  filter_acs_itt()  %>%
  group_by(CAssignment) %>%
  summarise(
    `Primary Outcome_No_`  = sum(1 - PO, na.rm = TRUE),
    `Primary Outcome_Yes_` = sum(PO, na.rm = TRUE),
    `Primary Outcome_Unknown_Total` = sum(is.na(PO)),
    `Primary Outcome_Unknown_Day 28 status` = sum(is.na(PO) & is.na(D28_death)),
    `Primary Outcome_Unknown_Vasopressor/inotropes` = 
      sum(is.na(PO) & !is.na(D28_death) & is.na(ANY_vasop))
  ) %>%
  gather(key, Frequency, -CAssignment, factor_key = T) %>%
  spread(CAssignment, Frequency) %>%
  separate(key, into = c("Measure", "Outcome", "Breakdown"), sep = "_") %>%
  mutate(
    Overall = C1 + C2 + C3 + C4,
    across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x[1:3])))
  )
kable(
  comp_table %>% select(-1),
  align = "llrr") %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Primary outcome", 1, nrow(comp_table))
Outcome Breakdown C1 C2 C3 C4 Overall
Primary outcome
No 561 (92) 576 (94) 259 (92) 43 (86) 1439 (92)
Yes 35 (6) 25 (4) 20 (7) 7 (14) 87 (6)
Unknown Total 14 (2) 12 (2) 4 (1) 0 (0) 30 (2)
Unknown Day 28 status 13 (2) 10 (2) 2 (1) 0 (0) 25 (2)
Unknown Vasopressor/inotropes 1 (0) 2 (0) 2 (1) 0 (0) 5 (0)
Table 1: Composite primary endpoint breakdown - overall.
Summarise the all components
tab <- bind_rows(comp_table, mort_table, vaso_table, vent_table)
tab_head <- all_dat %>% 
  filter_acs_itt() %>% 
  dplyr::count(CAssignment = factor(CAssignment, labels = intervention_labels()$CAssignment[-1])) %>%
  mutate(label = paste0(CAssignment, "<br>(n = ", n, ")")) %>%
  pull(label) %>%
  c(paste0("Overall<br><br>(n = ", nrow(all_dat %>% filter_acs_itt()), ")"))
colnames(tab)[4:8] <- tab_head
kable(
  tab %>% select(-1),
  align = "llrrrrr",
  escape = F) %>%
  kable_styling(font_size = 11) %>%
  collapse_rows(1, valign = "top", latex_hline = 'custom', custom_latex_hline = 1) %>%
  group_rows("Primary outcome", 1, nrow(comp_table)) %>%
  group_rows("Mortality", nrow(comp_table) + 1, nrow(comp_table) + nrow(mort_table)) %>%
  group_rows("Vasopressor/inotropes", 11, 15) %>%
  group_rows("Ventilation", 16, 20)
Outcome Breakdown Low
dose
(n = 610)
Intermediate
dose
(n = 613)
Low dose
with aspirin
(n = 283)
Therapeutic
dose
(n = 50)
Overall

(n = 1556)
Primary outcome
No 561 (92) 576 (94) 259 (92) 43 (86) 1439 (92)
Yes 35 (6) 25 (4) 20 (7) 7 (14) 87 (6)
Unknown Total 14 (2) 12 (2) 4 (1) 0 (0) 30 (2)
Unknown Day 28 status 13 (2) 10 (2) 2 (1) 0 (0) 25 (2)
Unknown Vasopressor/inotropes 1 (0) 2 (0) 2 (1) 0 (0) 5 (0)
Mortality
Alive at day 28 577 (95) 588 (96) 271 (96) 44 (88) 1480 (95)
Death within 28 days Total 19 (3) 15 (2) 10 (4) 6 (12) 50 (3)
Death within 28 days Prior to discharge 15 (2) 11 (2) 10 (4) 4 (8) 40 (3)
Death within 28 days Post-discharge 4 (1) 4 (1) 0 (0) 2 (4) 10 (1)
Unknown 14 (2) 10 (2) 2 (1) 0 (0) 26 (2)
Vasopressor/inotropes
Not required 591 (97) 595 (97) 273 (96) 48 (96) 1507 (97)
Use within 28 days Total 6 (1) 6 (1) 6 (2) 2 (4) 20 (1)
Use within 28 days Prior to discharge 5 (1) 6 (1) 6 (2) 2 (4) 19 (1)
Use within 28 days Post-discharge 2 (0) 0 (0) 0 (0) 1 (2) 3 (0)
Unknown 13 (2) 12 (2) 4 (1) 0 (0) 29 (2)
Ventilation
Not required 562 (92) 579 (94) 263 (93) 43 (86) 1447 (93)
Use within 28 days Total 35 (6) 24 (4) 18 (6) 7 (14) 84 (5)
Use within 28 days Prior to discharge 29 (5) 19 (3) 13 (5) 5 (10) 66 (4)
Use within 28 days Post-discharge 6 (1) 5 (1) 5 (2) 2 (4) 18 (1)
Unknown 13 (2) 10 (2) 2 (1) 0 (0) 25 (2)
Table 2: Composite primary endpoint breakdown.
Save table to outputs
colnames(tab)[4:8] <- linebreak(tab_head, align = "c", linebreaker = "<br>")
out <- kable(
  tab %>% select(-1),
  format = "latex",
  align = "llrrrrr",
  booktabs = TRUE,
  escape = F) %>%
  kable_styling(font_size = 9, latex_options = "HOLD_position") %>%
  collapse_rows(1, valign = "top", latex_hline = 'none') %>%
  group_rows("Primary outcome", 1, nrow(comp_table)) %>%
  group_rows("Mortality", nrow(comp_table) + 1, nrow(comp_table) + nrow(mort_table)) %>%
  group_rows("Vasopressor/inotropes", 11, 15) %>%
  group_rows("Ventilation", 16, 20)
save_tex_table(out, fn = "outcomes/primary/composite-breakdown")
Summarise the component combinations
levs <- c("None", "Death + Vas./Ino. + Ventilation", "Death + Vas./Ino.", "Death + Ventilation", "Vas./Ino. + Ventilation", 
          "Death", "Vas./Ino.", "Ventilation", "Unknown")
tab <- all_dat %>%
  filter_acs_itt()  %>%
  group_by(CAssignment) %>%
  summarise(
    `Composite outcomes` = case_when(
      PO == 0 ~ 0,
      is.na(PO) ~ 99,
      D28_death == 1 & ANY_vasop == 1 & ANY_vent == 1 ~ 1,
      D28_death == 1 & ANY_vasop == 1 & ANY_vent == 0 ~ 2,
      D28_death == 1 & ANY_vasop == 0 & ANY_vent == 1 ~ 3,
      D28_death == 0 & ANY_vasop == 1 & ANY_vent == 1 ~ 4,
      D28_death == 1 ~ 5,
      ANY_vasop == 1 ~ 6,
      ANY_vent == 1 ~ 7,
    )) %>%
  dplyr::count(CAssignment, `Composite outcomes`) %>%
  mutate(`Composite outcomes` = factor(`Composite outcomes`, levels = c(0:7, 99), labels = levs)) %>%
  spread(CAssignment, n, fill = 0) %>%
  mutate(Overall = C1 + C2 + C3 + C4) %>%
  mutate(across(C1:Overall, ~ sprintf("%i (%.0f)", .x, 100 * .x / sum(.x))))
colnames(tab)[-1] <- tab_head
colnames(tab)[1] <- "Composite outcomes, n (\\%)"
kable(
  tab ,
  align = "lrrrrr",
  escape = F) %>%
  kable_styling(font_size = 11)
Composite outcomes, n (\%) Low
dose
(n = 610)
Intermediate
dose
(n = 613)
Low dose
with aspirin
(n = 283)
Therapeutic
dose
(n = 50)
Overall

(n = 1556)
None 561 (92) 576 (94) 259 (92) 43 (86) 1439 (92)
Death + Vas./Ino. + Ventilation 6 (1) 3 (0) 3 (1) 2 (4) 14 (1)
Death + Ventilation 13 (2) 11 (2) 7 (2) 4 (8) 35 (2)
Vas./Ino. + Ventilation 0 (0) 2 (0) 1 (0) 0 (0) 3 (0)
Death 0 (0) 1 (0) 0 (0) 0 (0) 1 (0)
Vas./Ino. 0 (0) 1 (0) 2 (1) 0 (0) 3 (0)
Ventilation 16 (3) 7 (1) 7 (2) 1 (2) 31 (2)
Unknown 14 (2) 12 (2) 4 (1) 0 (0) 30 (2)
Table 3: Observed outcome combinations.
Code
colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
out <- kable(
  tab,
  format = "latex",
  align = "lrrrrr",
  booktabs = TRUE,
  linesep = "",
  escape = F) %>%
  kable_styling(
    font_size = 9, 
    latex_options = "HOLD_position")
save_tex_table(out, fn = "outcomes/primary/composite-combinations")

Descriptive Analyses

By Interim Analysis

Code
intdat <- all_dat %>%
  filter_acs_itt() %>%
  mutate(
    interim = as.character(case_when(
      RandDate < get_interim_dates()$meet_date[1] ~ 1,
      RandDate < get_interim_dates()$meet_date[2] ~ 2,
      RandDate < get_interim_dates()$meet_date[3] ~ 3,
      RandDate < get_interim_dates()$meet_date[4] ~ 4,
      TRUE ~ 4
    )),
    CAssignment = factor(
      CAssignment, labels = c("Low", "Intermediate", "Low\nwith aspirin", "Therapuetic"))
  )
intcounts <- intdat %>% 
  dplyr::count(interim) %>% 
  mutate(label = paste0("Interim ", interim, " (n = ", n, ")"))
trtcounts <- intdat %>% dplyr::count(CAssignment)
ovrtab <- intdat %>%
  group_by(interim = "Overall", CAssignment) %>%
  summarise(
    Randomised = 
      sprintf("%i", n()),
    Known = 
      sprintf("%i (%.0f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
    `Met primary outcome` = 
      sprintf("%i (%.0f)", sum(PO, na.rm = T), 100 * mean(PO, na.rm = T)),
    `Death` = 
      sprintf("%i (%.0f)", sum(D28_death, na.rm = T), 100 * mean(D28_death, na.rm = T)),
    `Vasopressor/inotropic support` = 
      sprintf("%i (%.0f)", sum(ANY_vasop, na.rm = T), 100 * mean(ANY_vasop, na.rm = T)),
    `Intensive respiratory support` = 
      sprintf("%i (%.0f)", sum(ANY_vent, na.rm = T), 100 * mean(ANY_vent, na.rm = T))
  ) %>%
  pivot_longer(`Randomised`:`Intensive respiratory support`, names_to = "Outcome") %>%
  pivot_wider(names_from = CAssignment, values_from = value, values_fill = "-")
inttab <- intdat %>%
  group_by(interim, CAssignment) %>%
  summarise(
    Randomised = 
      sprintf("%i", n()),
    Known = 
      sprintf("%i (%.0f)", sum(!is.na(PO)), 100 * mean(!is.na(PO))),
    `Met primary outcome` = 
      sprintf("%i (%.0f)", sum(PO, na.rm = T), 100 * mean(PO, na.rm = T)),
    `Death` = 
      sprintf("%i (%.0f)", sum(D28_death, na.rm = T), 100 * mean(D28_death, na.rm = T)),
    `Vasopressor/inotropic support` = 
      sprintf("%i (%.0f)", sum(ANY_vasop, na.rm = T), 100 * mean(ANY_vasop, na.rm = T)),
    `Intensive respiratory support` = 
      sprintf("%i (%.0f)", sum(ANY_vent, na.rm = T), 100 * mean(ANY_vent, na.rm = T))
  ) %>%
  pivot_longer(`Randomised`:`Intensive respiratory support`, names_to = "Outcome") %>%
  pivot_wider(names_from = CAssignment, values_from = value, values_fill = "-")
tab <- bind_rows(ovrtab, inttab)

savetab <- kable(
  tab[, -1],
  format = "latex",
  align = "lrrrr",
  booktabs = TRUE
) %>%
  kable_styling(
    font_size = 9,
    latex_options = "HOLD_position"
  ) %>%
  group_rows(paste0("Overall (n = ", sum(intcounts$n), ")"), 1, 6) %>%
  group_rows(intcounts$label[1], 7, 12) %>%
  group_rows(intcounts$label[2], 13, 18) %>%
  group_rows(intcounts$label[3], 19, 24) %>%
  group_rows(intcounts$label[4], 25, 30)
save_tex_table(savetab, fn = "outcomes/primary/primary-composite-interims")

kable(
  tab[, -1],
  format = "html",
  align = "lrrrr",
) %>%
  kable_styling(
    font_size = 9,
    latex_options = "HOLD_position"
  ) %>%
  group_rows(paste0("Overall (n = ", sum(intcounts$n), ")"), 1, 6) %>%
  group_rows(intcounts$label[1], 7, 12) %>%
  group_rows(intcounts$label[2], 13, 18) %>%
  group_rows(intcounts$label[3], 19, 24) %>%
  group_rows(intcounts$label[4], 25, 30)
Outcome Low Intermediate Low with aspirin Therapuetic
Overall (n = 1556)
Randomised 610 613 283 50
Known 596 (98) 601 (98) 279 (99) 50 (100)
Met primary outcome 35 (6) 25 (4) 20 (7) 7 (14)
Death 19 (3) 15 (2) 10 (4) 6 (12)
Vasopressor/inotropic support 6 (1) 6 (1) 6 (2) 2 (4)
Intensive respiratory support 35 (6) 24 (4) 18 (6) 7 (14)
Interim 1 (n = 685)
Randomised 234 226 225 -
Known 227 (97) 222 (98) 222 (99) -
Met primary outcome 22 (10) 15 (7) 19 (9) -
Death 13 (6) 11 (5) 9 (4) -
Vasopressor/inotropic support 2 (1) 3 (1) 6 (3) -
Intensive respiratory support 22 (10) 15 (7) 17 (8) -
Interim 2 (n = 246)
Randomised 90 98 58 -
Known 90 (100) 97 (99) 57 (98) -
Met primary outcome 2 (2) 3 (3) 1 (2) -
Death 1 (1) 3 (3) 1 (2) -
Vasopressor/inotropic support 1 (1) 0 (0) 0 (0) -
Intensive respiratory support 2 (2) 3 (3) 1 (2) -
Interim 3 (n = 384)
Randomised 177 172 - 35
Known 173 (98) 167 (97) - 35 (100)
Met primary outcome 8 (5) 4 (2) - 6 (17)
Death 4 (2) 0 (0) - 5 (14)
Vasopressor/inotropic support 3 (2) 2 (1) - 2 (6)
Intensive respiratory support 8 (5) 4 (2) - 6 (17)
Interim 4 (n = 241)
Randomised 109 117 - 15
Known 106 (97) 115 (98) - 15 (100)
Met primary outcome 3 (3) 3 (3) - 1 (7)
Death 1 (1) 1 (1) - 1 (7)
Vasopressor/inotropic support 0 (0) 1 (1) - 0 (0)
Intensive respiratory support 3 (3) 2 (2) - 1 (7)

Age

Relationship age to outcome
agedat <- acs_itt_dat %>%
  dplyr::count(PO, AgeAtEntry) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nwith outcome", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-age.pdf"), p, height = 2, width = 6)
p

Relationship age 265 60 to outcome
acs_itt_dat %>%
  dplyr::count(`Age 60+` = agegte60, PO) %>%
  group_by(`Age 60+`) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  kable(digits = 2) %>%
  kable_styling(bootstrap_options = "striped")
Age 60+ 0 1 <NA> p_1 p_miss
0 1034 47 21 0.04 0.02
1 405 40 9 0.09 0.02
Table 4: Relationship between age >= 60 and the primary outcome.

Country

Relationship country to outcome
tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")), PO) %>%
  group_by(Country) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion meeting\nprimary outcome", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-country.pdf"), p, height = 2, width = 6)
p

Site

Relationship site to outcome
tdat <- all_dat %>%
  filter_acs_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    PO) %>%
  group_by(Country, Site) %>%
  spread(PO, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()

p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion meeting\nprimary outcome", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-country-site.pdf"), p, height = 4, width = 6.25)
p

Calendar Time

Relationship calendar date to outcome
caldat <- acs_itt_dat %>%
  dplyr::count(PO, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion meeting\nprimary outcome", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, 0.25)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-calendar-time.pdf"), p, height = 2, width = 6)
p

Days since symptom onset

Relationship days since sympton onset to outcome
dsfsdat <- acs_itt_dat %>%
  dplyr::count(PO, dsfs) %>% 
  spread(PO, n, fill = 0) %>% 
  mutate(n = `1` + `0` + `<NA>`,
         p = `1` / (`1` + `0`))
dsfsmod <- glm(
  cbind(`1`, `0`) ~ dsfs, 
  data = dsfsdat, 
  family = binomial())
dsfsdat <- dsfsdat %>%
  mutate(
    ypred = predict(dsfsmod, newdata = dsfsdat, type = "response")
  )
p1 <- ggplot(dsfsdat, aes(dsfs, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 7.5, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Days since first symptoms")
p2 <- ggplot(dsfsdat, aes(dsfs, p)) +
    geom_point() +
    geom_vline(xintercept = 7.5, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\nwith outcome", 
         x = "Days since first symptoms")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "primary")
ggsave(file.path(path, "outcome-dsfs.pdf"), p, height = 2, width = 6)
p

Anti-coagulation Domain

Relationship anti-coagulation to outcome
make_po_table(
    all_dat %>%
      filter_acs_itt()
)
n (\%) Low
dose
Intermediate
dose
Low dose
with aspirin
Therapeutic
dose
Overall
Randomised 610 613 283 50 1556
Outcome missing 14 (2.3) 12 (2.0) 4 (1.4) 0 (0.0) 30 (1.9)
Outcome observed 596 (97.7) 601 (98.0) 279 (98.6) 50 (100.0) 1526 (98.1)
Met primary outcome 35 (5.9) 25 (4.2) 20 (7.2) 7 (14.0) 87 (5.7)
Table 5: Primary outcome by anti-coagulation intervention.
Relationship anti-coagulation to outcome
save_tex_table(make_po_table(
      all_dat %>%
      filter_acs_itt(),
      "latex"), 
      "outcomes/primary/anticoagulation-summary")

Analyses

ACS-ITT

In these analyses, the population is restricted to participants who were randomised to the anti-coagulation domain.

Treatment only

Data and sampling treatment only model
res <- fit_primary_model(acs_itt_nona_dat, logistic_mod, seed = 17250)
names(res$drws$OR) <- intervention_labels2()$CAssignment[-(1:2)]
Code
plot_or_densities(res$drws$OR)

Code
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.71 (0.43, 1.20) 0.74 (0.20) 0.90
Low-dose with aspirin 1.23 (0.70, 2.14) 1.28 (0.37) 0.23
Therapeutic-dose 2.40 (0.95, 5.34) 2.60 (1.16) 0.03
Table 6: Posterior summaries for odds ratio contrasts.
Code
res$fit$summary("beta")
# A tibble: 6 × 10
  variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
  <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
1 beta[1]  -2.62   -2.61   0.143 0.143 -2.86   -2.39    1.00   15499.   14290.
2 beta[2]   0.485   0.494  0.334 0.330 -0.0778  1.02    1.00   17188.   14447.
3 beta[3]  -0.210  -0.210  0.210 0.211 -0.554   0.136   1.00   17721.   15370.
4 beta[4]  -0.695  -0.694  0.246 0.245 -1.10   -0.289   1.00   17656.   15552.
5 beta[5]  -0.0510 -0.0333 0.399 0.396 -0.735   0.580   1.00   20107.   14608.
6 beta[6]  -0.785  -0.774  0.504 0.501 -1.63    0.0235  1.00   21997.   13319.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0843 1.0378 1.0546 1.1228 1.0945 1.0608 0.9765 1.0343
Code
mcmc_trace(res$drws["beta"])

Treatment plus covariates (no site or epoch)

Under this model a component is added for country (region), \(\mathbb{I}[\text{age}\geq 60]\), and intervention specific ineligibility (C3 only)

New terms:

\[ \begin{aligned} \rho_1 &= 0 \\ \rho_r &\overset{\text{iid}}{\sim} \text{Normal}(0, 1), r=2,...,R,\\ \beta_{\text{age}\geq 60} &\sim \text{Normal}(0, 2.5^2) \\ \beta_{\text{inelgC3}} &\sim \text{Normal}(0, 10) \end{aligned} \]

Code
res <- fit_primary_model(
  acs_itt_nona_dat, 
  logistic_mod, 
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 71236)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")
Code
plot_or_densities(res$drws$OR)

Code
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.71 (0.42, 1.21) 0.74 (0.20) 0.90
Low-dose with aspirin 1.32 (0.73, 2.33) 1.37 (0.41) 0.18
Therapeutic-dose 1.85 (0.68, 4.60) 2.05 (1.03) 0.11
Ineligible aspirin 2.39 (0.74, 6.29) 2.67 (1.44) 0.07
Age ≥ 60 1.96 (1.25, 3.05) 2.01 (0.46) 0.00
Australia/New Zealand 1.00 (0.29, 3.06) 1.17 (0.74) 0.50
Nepal 1.43 (0.65, 2.92) 1.53 (0.59) 0.18
Table 7: Posterior summaries for odds ratio contrasts.
Code
res$fit$summary("beta")
# A tibble: 10 × 10
   variable    mean   median    sd   mad      q5      q95  rhat ess_bulk
   <chr>      <dbl>    <dbl> <dbl> <dbl>   <dbl>    <dbl> <dbl>    <dbl>
 1 beta[1]  -3.00   -3.00    0.198 0.197 -3.33   -2.68     1.00   14028.
 2 beta[2]   0.253   0.255   0.378 0.378 -0.379   0.866    1.00   15820.
 3 beta[3]  -0.155  -0.158   0.216 0.216 -0.506   0.198    1.00   19065.
 4 beta[4]  -0.626  -0.628   0.261 0.262 -1.05   -0.193    1.00   17111.
 5 beta[5]   0.0869  0.0802  0.646 0.645 -0.960   1.15     1.00   16336.
 6 beta[6]  -0.795  -0.787   0.497 0.495 -1.62    0.00501  1.00   21338.
 7 beta[7]   0.843   0.873   0.542 0.532 -0.0981  1.68     1.00   22546.
 8 beta[8]   0.674   0.675   0.226 0.224  0.298   1.05     1.00   18977.
 9 beta[9]  -0.0176 -0.00389 0.606 0.600 -1.03    0.953    1.00   16486.
10 beta[10]  0.352   0.358   0.382 0.381 -0.290   0.965    1.00   18103.
# … with 1 more variable: ess_tail <dbl>
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0215 1.0098 1.0926 1.0081 1.0257 0.9988 1.0248 1.0574
Code
mcmc_trace(res$drws["beta"])

Treatment plus covariates plus site (no epoch)

Under this model a component is added for site nested within country.

New terms

\[ \begin{aligned} \xi_{rs} &\overset{\text{iid}}{\sim}\text{Normal}(0,\sigma^2_{r}), s=1,...,S_r, r=1,...,R \\ \sigma_r &\overset{\text{iid}}{\sim} \text{Half-}t(3, 1). \end{aligned} \]

Data and sampling no site or epoch model
res <- fit_primary_model(
  acs_itt_nona_dat, 
  logistic_site_mod, 
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 52765,
  adapt_delta = 0.95)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")
Code
plot_or_densities(res$drws$OR)

Code
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.72 (0.41, 1.24) 0.75 (0.21) 0.88
Low-dose with aspirin 0.94 (0.51, 1.74) 0.99 (0.31) 0.57
Therapeutic-dose 2.07 (0.73, 5.66) 2.37 (1.31) 0.08
Ineligible aspirin 2.72 (0.75, 8.52) 3.20 (2.07) 0.06
Age ≥ 60 1.80 (1.11, 2.93) 1.86 (0.46) 0.01
Australia/New Zealand 1.01 (0.22, 4.21) 1.32 (1.10) 0.49
Nepal 1.87 (0.46, 6.63) 2.30 (1.70) 0.17
Table 8: Posterior summaries for odds ratio contrasts.
Code
plot_site_terms(
  exp(res$drws$gamma_site), 
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)

Code
res$fit$summary(c("beta", "gamma_site", "tau_site"))
# A tibble: 39 × 10
   variable     mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>       <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -3.56    -3.55   0.532 0.519 -4.45   -2.71    1.00    5496.    8830.
 2 beta[2]   0.576    0.577  0.404 0.399 -0.101   1.23    1.00   17006.   14724.
 3 beta[3]  -0.0975  -0.0956 0.228 0.228 -0.474   0.280   1.00   18674.   15698.
 4 beta[4]  -0.526   -0.526  0.272 0.272 -0.974  -0.0808  1.00   17773.   14426.
 5 beta[5]   0.195    0.196  0.700 0.702 -0.961   1.35    1.00   16795.   14547.
 6 beta[6]  -0.966   -0.955  0.547 0.542 -1.88   -0.0897  1.00   23767.   14794.
 7 beta[7]   0.978    1.00   0.619 0.613 -0.0588  1.95    1.00   22493.   14369.
 8 beta[8]   0.588    0.586  0.246 0.243  0.184   0.992   1.00   25540.   15554.
 9 beta[9]   0.00242  0.0117 0.750 0.745 -1.26    1.20    1.00   12505.   14538.
10 beta[10]  0.614    0.627  0.667 0.645 -0.510   1.68    1.00    9105.   12175.
# … with 29 more rows
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 2

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8537 0.8557 0.8084 0.8207 0.7391 0.8056 0.8420 0.7953
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["tau_site"])

Treatment plus covariates plus epoch (no site)

Replace site terms with epoch terms

\[ \begin{aligned} \tau_1 &= 0 \\ \tau_t &= \tau_{t-1} + \sigma_\tau\epsilon_t \\ \epsilon_t &\overset{\text{iid}}{\sim} \text{Normal}(0,1) \\ \sigma_\tau &\sim \text{Half-}t(3, 1) \end{aligned} \]

Data and sampling no site model
res <- fit_primary_model(
  acs_itt_nona_dat, 
  logistic_epoch_mod, 
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 76405,
  adapt_delta = 0.98)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")
Code
plot_or_densities(res$drws$OR)

Code
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.73 (0.43, 1.24) 0.76 (0.21) 0.87
Low-dose with aspirin 0.92 (0.50, 1.68) 0.97 (0.31) 0.60
Therapeutic-dose 2.30 (0.83, 6.08) 2.59 (1.39) 0.05
Ineligible aspirin 1.70 (0.52, 4.59) 1.92 (1.07) 0.17
Age ≥ 60 1.95 (1.22, 3.09) 2.00 (0.48) 0.00
Australia/New Zealand 1.41 (0.40, 4.58) 1.70 (1.11) 0.29
Nepal 2.91 (1.15, 7.19) 3.23 (1.58) 0.01
Table 9: Posterior summaries for odds ratio contrasts.
Code
plot_epoch_terms(
  exp(res$drws$gamma_epoch)
)

Code
res$fit$summary(c("beta", "gamma_epoch", "tau_epoch"))
# A tibble: 23 × 10
   variable   mean median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>     <dbl>  <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.22  -4.19  0.528 0.511 -5.13   -3.39    1.00   15486.   12840.
 2 beta[2]   0.659  0.662 0.406 0.402 -0.0127  1.32    1.00   18941.   15955.
 3 beta[3]  -0.125 -0.122 0.219 0.219 -0.486   0.235   1.00   20213.   15097.
 4 beta[4]  -0.535 -0.534 0.266 0.267 -0.976  -0.102   1.00   18884.   15031.
 5 beta[5]   0.348  0.340 0.654 0.651 -0.709   1.43    1.00   16689.   14620.
 6 beta[6]  -0.833 -0.819 0.512 0.507 -1.70   -0.0150  1.00   25915.   14359.
 7 beta[7]   0.504  0.530 0.553 0.547 -0.431   1.37    1.00   24106.   13640.
 8 beta[8]   0.663  0.667 0.237 0.236  0.273   1.05    1.00   25646.   14940.
 9 beta[9]   0.339  0.346 0.622 0.625 -0.707   1.34    1.00   17486.   14074.
10 beta[10]  1.06   1.07  0.465 0.464  0.291   1.82    1.00   17849.   14963.
# … with 13 more rows
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8896 0.7985 0.8545 0.8386 0.8910 0.8580 0.8231 0.8728
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Code
mcmc_trace(res$drws["tau_epoch"])

Primary Model

Data and sampling primary model
res <- fit_primary_model(
  acs_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
res$fit$save_object(file.path("outputs", "models", "primary", "acs-itt-primary.rds"))
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")
Code
save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-summary-table")
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.74 (0.43, 1.27) 0.77 (0.22) 0.86
Low-dose with aspirin 0.88 (0.47, 1.64) 0.93 (0.30) 0.65
Therapeutic-dose 2.22 (0.77, 6.20) 2.55 (1.45) 0.07
Ineligible aspirin 2.50 (0.69, 7.93) 2.96 (1.94) 0.08
Age ≥ 60 1.88 (1.16, 3.04) 1.94 (0.49) 0.01
Australia/New Zealand 1.05 (0.23, 4.39) 1.37 (1.13) 0.47
Nepal 2.17 (0.51, 7.67) 2.66 (1.97) 0.13
Posterior summaries for model parameters.
Code
save_tex_table(
  decision_quantity_summary_table(c("Low-dose" = rvar(1), res$drws$OR[1:3]), "latex"),
  "outcomes/primary/primary-model-acs-itt-decision-table")
decision_quantity_summary_table(
  c("Low-dose" = rvar(1), res$drws$OR[1:3])
)
Intervention Posterior Superior Pr(OR = min(OR)) Effective Pr(OR < 1) Futile Pr(OR > 1/1.1) Equivalent Pr(1/1.1 < OR < 1.1)
Low-dose 1.00 0.08 - - -
Intermediate-dose 0.74 (0.43, 1.27) 0.63 0.86 0.23 0.15
Low-dose with aspirin 0.88 (0.47, 1.64) 0.27 0.65 0.47 0.22
Therapeutic-dose 2.22 (0.77, 6.20) 0.02 0.07 0.95 0.04
Table 10: Decision quantity summaries
Code
p <- plot_or_densities(res$drws$OR[1:3]) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "primary", "primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Code
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "primary", "primary-model-acs-itt-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Code
res$fit$summary(c("beta", "gamma_epoch", "tau_epoch", "gamma_site", "tau_site"))
# A tibble: 52 × 10
   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.07   -4.06   0.671 0.654 -5.21   -3.01    1.00    8389.   12597.
 2 beta[2]   0.666   0.668  0.420 0.416 -0.0262  1.36    1.00   20492.   15814.
 3 beta[3]  -0.105  -0.103  0.228 0.227 -0.483   0.270   1.00   22714.   15494.
 4 beta[4]  -0.502  -0.499  0.277 0.277 -0.961  -0.0501  1.00   20392.   14539.
 5 beta[5]   0.224   0.218  0.698 0.697 -0.910   1.37    1.00   19342.   15216.
 6 beta[6]  -0.947  -0.944  0.542 0.533 -1.84   -0.0635  1.00   27617.   15229.
 7 beta[7]   0.899   0.915  0.624 0.618 -0.162   1.89    1.00   27309.   14051.
 8 beta[8]   0.634   0.633  0.248 0.249  0.222   1.04    1.00   27924.   14818.
 9 beta[9]   0.0392  0.0515 0.755 0.752 -1.24    1.26    1.00   14945.   14640.
10 beta[10]  0.749   0.773  0.687 0.683 -0.428   1.83    1.00   11783.   13803.
# … with 42 more rows
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8537 0.7930 0.7605 0.7244 0.7954 0.7820 0.7752 0.7430
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Code
mcmc_trace(res$drws["tau_epoch"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["tau_site"])

Posterior Predictive Checks

Code
y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc))

grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(PO),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}

ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = country, grp = site_raw) %>%
  summarise(
    mean_y = mean(PO),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )

p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "primary",
                 "primary-model-acs-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Sensitivity Analyses

Sens: Treatment Contrasts

Code
res <- fit_primary_model(
  acs_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  ctr = contr.treatment,
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")
Code
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.71 (0.42, 1.21) 0.74 (0.20) 0.90
Low-dose with aspirin 0.85 (0.46, 1.55) 0.89 (0.28) 0.70
Therapeutic-dose 2.04 (0.72, 5.39) 2.29 (1.22) 0.08
Ineligible aspirin 2.44 (0.66, 7.74) 2.89 (1.92) 0.08
Age ≥ 60 1.90 (1.17, 3.07) 1.96 (0.49) 0.01
Australia/New Zealand 1.18 (0.25, 4.84) 1.52 (1.27) 0.41
Nepal 2.10 (0.50, 7.44) 2.58 (1.93) 0.14
Posterior summaries for model parameters under treatment contrast priors.
Code
save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-trt-ctr-summary-table")
Code
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "primary",
                 "primary-model-acs-itt-trt-ctr-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Weakly Informative Priors

As a prior sensitivity assessment, the primary model is re-fit under weakly informative priors (Normal(0, 10)) for all fixed-effect parameters.

Code
res <- fit_primary_model(
  acs_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 10, 10, 10),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")
Code
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.74 (0.43, 1.29) 0.77 (0.22) 0.86
Low-dose with aspirin 0.89 (0.47, 1.65) 0.93 (0.30) 0.64
Therapeutic-dose 2.10 (0.73, 5.83) 2.41 (1.36) 0.08
Ineligible aspirin 2.51 (0.67, 8.19) 2.99 (2.00) 0.08
Age ≥ 60 1.89 (1.16, 3.10) 1.95 (0.50) 0.01
Australia/New Zealand 1.37 (0.09, 12.45) 2.59 (4.28) 0.40
Nepal 4.05 (0.54, 32.02) 7.29 (17.05) 0.07
Posterior summaries for model parameters under treatment contrast priors.
Code
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
p

Sens: FAS-ITT

In these analyses, the population is includes all participants who were randomised to the platform, even if not randomised into the anticoagulation domain.

Code
res <- fit_primary_model(
  fas_itt_nona_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.74 (0.43, 1.27) 0.77 (0.22) 0.86
Low-dose with aspirin 0.89 (0.48, 1.62) 0.93 (0.29) 0.64
Therapeutic-dose 2.17 (0.77, 6.00) 2.48 (1.37) 0.07
Ineligible aspirin 2.16 (0.60, 6.75) 2.55 (1.64) 0.11
Age ≥ 60 1.89 (1.18, 3.05) 1.95 (0.48) 0.01
Australia/New Zealand 1.11 (0.26, 4.56) 1.43 (1.16) 0.44
Nepal 2.19 (0.50, 7.89) 2.70 (2.00) 0.13
Table 11: Posterior summaries for odds ratio contrasts.
Code
save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-fas-itt-summary-table")
Code
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "primary", "primary-model-fas-itt-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Worst-Case

Under this scenario, all participants with unknown primary outcome status were were to assumed to have met the primary outcome.

Code
make_po_table(acs_itt_wc_dat)
n (\%) Low
dose
Intermediate
dose
Low dose
with aspirin
Therapeutic
dose
Overall
Randomised 610 613 283 50 1556
Outcome missing 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0)
Outcome observed 610 (100.0) 613 (100.0) 283 (100.0) 50 (100.0) 1556 (100.0)
Met primary outcome 49 (8.0) 37 (6.0) 24 (8.5) 7 (14.0) 117 (7.5)
Code
save_tex_table(
  make_po_table(acs_itt_wc_dat, "latex"), 
  "outcomes/primary/acs-itt-worst-case-anticoagulation-summary")
Code
res <- fit_primary_model(
  acs_itt_wc_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code
save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-worst-case-summary-table")
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.74 (0.46, 1.19) 0.77 (0.19) 0.89
Low-dose with aspirin 0.81 (0.46, 1.40) 0.84 (0.24) 0.77
Therapeutic-dose 1.37 (0.50, 3.51) 1.53 (0.78) 0.26
Ineligible aspirin 1.22 (0.35, 3.67) 1.42 (0.89) 0.37
Age ≥ 60 1.67 (1.07, 2.59) 1.71 (0.39) 0.01
Australia/New Zealand 2.37 (0.54, 9.18) 2.99 (2.34) 0.12
Nepal 1.89 (0.43, 6.83) 2.31 (1.72) 0.18
Posterior summaries for model parameters.
Code
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal")))
pth <- file.path("outputs", "figures", "outcomes", "primary", 
                 "primary-model-acs-itt-worst-case-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Best-Case

Under this scenario, all participants with unknown primary outcome status were were to assumed to have not met the primary outcome.

Code
save_tex_table(
  make_po_table(acs_itt_bc_dat, "latex"), 
  "outcomes/primary/acs-itt-best-case-anticoagulation-summary")
make_po_table(acs_itt_bc_dat)
n (\%) Low
dose
Intermediate
dose
Low dose
with aspirin
Therapeutic
dose
Overall
Randomised 610 613 283 50 1556
Outcome missing 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0) 0 (0.0)
Outcome observed 610 (100.0) 613 (100.0) 283 (100.0) 50 (100.0) 1556 (100.0)
Met primary outcome 35 (5.7) 25 (4.1) 20 (7.1) 7 (14.0) 87 (5.6)
Code
res <- fit_primary_model(
  acs_itt_bc_dat,
  primary_mod,
  vars = c("inelgc3", "agegte60", "ctry"),
  beta_sd_var = c(10, 2.5, 1, 1),
  seed = 32915,
  adapt_delta = 0.99
)
names(res$drws$OR) <- c(
  intervention_labels2()$CAssignment[-(1:2)], 
  "Ineligible aspirin", "Age \u2265 60", "Australia/New Zealand", "Nepal")

Code
save_tex_table(
  odds_ratio_summary_table(res$drws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-best-case-summary-table")
odds_ratio_summary_table(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.75 (0.43, 1.30) 0.78 (0.22) 0.85
Low-dose with aspirin 0.91 (0.49, 1.67) 0.95 (0.31) 0.62
Therapeutic-dose 2.32 (0.80, 6.45) 2.66 (1.49) 0.06
Ineligible aspirin 2.83 (0.80, 9.03) 3.37 (2.19) 0.05
Age ≥ 60 1.86 (1.13, 3.04) 1.92 (0.49) 0.01
Australia/New Zealand 1.00 (0.22, 4.07) 1.29 (1.06) 0.50
Nepal 2.15 (0.49, 7.71) 2.64 (1.95) 0.14
Posterior summaries for model parameters.
Code
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch), 
  exp(res$drws$gamma_site),
  factor(
    res$dat$region_by_site, 
    labels = c("India", "Australia\nNew Zealand", "Nepal")))
pth <- file.path("outputs", "figures", "outcomes", "primary", 
                 "primary-model-acs-itt-best-case-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Sens: Concurrent Randomisation

Intermediate-dose

Code
save_tex_table(
  make_po_table(
    acs_itt_concurc2_dat,
    "latex"), 
  "outcomes/primary/acs-itt-c2-summary")
make_po_table(acs_itt_concurc2_dat)
n (\%) Low
dose
Intermediate
dose
Overall
Randomised 610 613 1223
Outcome missing 14 (2.3) 12 (2.0) 26 (2.1)
Outcome observed 596 (97.7) 601 (98.0) 1197 (97.9)
Met primary outcome 35 (5.9) 25 (4.2) 60 (5.0)
Table 12: Primary outcome summary for participants randomised to either low-dose or intermediate dose arms (ACS-ITT-intermediate).
Code
ctr <- contr.orthonorm
X <- model.matrix(
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc2_nona_dat,
  contrasts.arg = list(CAssignment = ctr))
y <- acs_itt_concurc2_nona_dat$PO
sdat <- list(
  N = nrow(X),
  K = ncol(X),
  X = X,
  y = y,
  beta_sd = c(2.5, 1, 2.5, 1, 1)
)
snk <- capture.output(
  sfit <- logistic_mod$sample(
    data = sdat,
    seed = 3274,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500
  )
)
sdrws <- as_draws_rvars(sfit$draws("beta"))
names(sdrws$beta) <- colnames(X)
sdrws$beta_trt <- as.vector(
  transform_coding(cbind(1, ctr(2)), cbind(1, contr.treatment(2))) %**% sdrws$beta[1:2]
)
sdrws$OR <- c(
  "Intermediate-dose" = exp(sdrws$beta_trt[2]),
  "Age \u2265 60" = exp(sdrws$beta[3]),
  "Australia/New Zealand" = exp(sdrws$beta[4]),
  "Nepal" = exp(sdrws$beta[5])
)
Code
save_tex_table(
  odds_ratio_summary_table(sdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-c2-summary-table")
odds_ratio_summary_table(sdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.69 (0.41, 1.16) 0.71 (0.19) 0.92
Age ≥ 60 2.19 (1.29, 3.66) 2.26 (0.61) 0.00
Australia/New Zealand 1.39 (0.64, 2.77) 1.47 (0.55) 0.20
Nepal 1.06 (0.40, 2.38) 1.15 (0.51) 0.45
Table 13: Parameter posterior summaries.
Code
res$fit$summary(c("beta"))
# A tibble: 10 × 10
   variable    mean   median    sd   mad      q5     q95  rhat ess_bulk ess_tail
   <chr>      <dbl>    <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
 1 beta[1]  -4.08   -4.07    0.667 0.652 -5.20   -3.01    1.00    8163.   11086.
 2 beta[2]   0.678   0.683   0.420 0.417 -0.0206  1.37    1.00   22134.   15754.
 3 beta[3]  -0.130  -0.131   0.229 0.227 -0.508   0.245   1.00   22749.   15590.
 4 beta[4]  -0.517  -0.517   0.275 0.276 -0.969  -0.0673  1.00   21721.   16386.
 5 beta[5]   0.258   0.261   0.698 0.698 -0.898   1.41    1.00   21338.   14916.
 6 beta[6]  -0.979  -0.968   0.542 0.542 -1.89   -0.112   1.00   30272.   14314.
 7 beta[7]   1.03    1.04    0.619 0.616 -0.0154  2.02    1.00   29722.   14530.
 8 beta[8]   0.620   0.620   0.252 0.252  0.204   1.03    1.00   32053.   14952.
 9 beta[9]  -0.0154 -0.00183 0.743 0.748 -1.26    1.18    1.00   14921.   15151.
10 beta[10]  0.742   0.767   0.692 0.669 -0.432   1.83    1.00   11785.   14311.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7842 0.7528 0.7704 0.8127 0.7986 0.7521 0.7614 0.7696
Code
mcmc_trace(res$drws["beta"])

Low-dose with aspirin

Code
save_tex_table(
  make_po_table(
    acs_itt_concurc3_dat,
    "latex"), 
  "outcomes/primary/acs-itt-c3-summary")
make_po_table(acs_itt_concurc3_dat)
n (\%) Low
dose
Intermediate
dose
Low dose
with aspirin
Overall
Randomised 299 298 283 880
Outcome missing 7 (2.3) 5 (1.7) 4 (1.4) 16 (1.8)
Outcome observed 292 (97.7) 293 (98.3) 279 (98.6) 864 (98.2)
Met primary outcome 21 (7.2) 16 (5.5) 20 (7.2) 57 (6.6)
Table 14: Primary outcome summary for participants randomised whilst the low-dose plus aspirin arm was available (ACS-ITT-aspirin).
Code
X <- model.matrix(
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc3_nona_dat,
  contrasts.arg = list(CAssignment = contr.orthonorm))
y <- acs_itt_concurc3_nona_dat$PO
sdat <- list(
  N = nrow(X),
  K = ncol(X),
  X = X,
  y = y,
  beta_sd = c(2.5, 1, 1, 2.5, 1)
)
snk <- capture.output(
  sfit <- logistic_mod$sample(
    data = sdat,
    seed = 3274,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500
  )
)
sdrws <- as_draws_rvars(sfit$draws("beta"))
names(sdrws$beta) <- colnames(X)

sdrws$beta_trt <- as.vector(
  transform_coding(cbind(1, contr.orthonorm(3)), cbind(1, contr.treatment(3))) %**% sdrws$beta[1:3]
)

sdrws$OR <- c(
  "Intermediate-dose" = exp(sdrws$beta_trt[2]),
  "Low-dose with aspirin" = exp(sdrws$beta_trt[3]),
  "Age \u2265 60" = exp(sdrws$beta[4]),
  "Australia/New Zealand" = exp(sdrws$beta[5])
)
Code
save_tex_table(
  odds_ratio_summary_table(sdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-c3-summary-table")
odds_ratio_summary_table(sdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate-dose 0.75 (0.39, 1.45) 0.79 (0.28) 0.80
Low-dose with aspirin 0.93 (0.50, 1.73) 0.98 (0.32) 0.59
Age ≥ 60 1.91 (1.09, 3.28) 1.98 (0.57) 0.01
Australia/New Zealand 0.31 (0.07, 1.09) 0.38 (0.27) 0.97
Table 15: Parameter posterior summaries.

Therapeutic-dose

Code
save_tex_table(
  make_po_table(
    acs_itt_concurc4_dat,
    "latex"), 
  "outcomes/primary/acs-itt-c4-summary")
make_po_table(acs_itt_concurc4_dat)
n (\%) Low
dose
Intermediate
dose
Therapeutic
dose
Overall
Randomised 79 65 50 194
Outcome missing 4 (5.1) 3 (4.6) 0 (0.0) 7 (3.6)
Outcome observed 75 (94.9) 62 (95.4) 50 (100.0) 187 (96.4)
Met primary outcome 6 (8.0) 3 (4.8) 7 (14.0) 16 (8.6)
Table 16: Primary outcome summary for participants randomised under protocol 5.0 (ACS-ITT-therapeutic).
Code
X <- model.matrix(
  ~ CAssignment + agegte60 + ctry,
  data = acs_itt_concurc4_nona_dat,
  contrasts.arg = list(CAssignment = contr.orthonorm))
y <- acs_itt_concurc4_nona_dat$PO
sdat <- list(
  N = nrow(X),
  K = ncol(X),
  X = X,
  y = y,
  beta_sd = c(2.5, 1, 1, 2.5, 1, 1)
)
snk <- capture.output(
  sfit <- logistic_mod$sample(
    data = sdat,
    seed = 3274,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500
  )
)
sdrws <- as_draws_rvars(sfit$draws("beta"))
names(sdrws$beta) <- colnames(X)

sdrws$beta_trt <- as.vector(
  transform_coding(cbind(1, contr.orthonorm(3)), cbind(1, contr.treatment(3))) %**% sdrws$beta[1:3]
)

sdrws$OR <- c(
  "Intermediate dose" = exp(sdrws$beta_trt[2]),
  "Therapeutic dose" = exp(sdrws$beta_trt[3]),
  "Age \u2265 60" = exp(sdrws$beta[4]),
  "India" = exp(sdrws$beta[5]),
  "Australia/New Zealand" = exp(sdrws$beta[6])
)
Code
save_tex_table(
  odds_ratio_summary_table(sdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-itt-c4-summary-table")
odds_ratio_summary_table(sdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate dose 0.65 (0.18, 2.25) 0.80 (0.57) 0.75
Therapeutic dose 1.63 (0.54, 5.07) 1.93 (1.20) 0.19
Age ≥ 60 1.88 (0.66, 5.45) 2.18 (1.27) 0.12
India 0.39 (0.08, 1.62) 0.51 (0.41) 0.90
Australia/New Zealand 1.05 (0.37, 2.84) 1.19 (0.65) 0.46
Table 17: Parameter posterior summaries.

Per-Protocol Analysis

In the per-protocol analysis, participants who did not satisfy the intervention protocol (as determined by blinded review) are excluded from the analysis set.

Code
all_pp %>%
  filter(ENR_rec == 1) %>%
  dplyr::count(Reason)
# A tibble: 4 × 2
  Reason                              n
  <chr>                           <int>
1 Patient Withdrawal from Domain     26
2 Protocol Deviation                 28
3 Protocol Deviation & Withdrawal     1
4 <NA>                             1544
Code
acs_itt_dat %>%
  left_join(pp, by = "StudyPatientID") %>%
  dplyr::count(CAssignment, Reason) %>%
  spread(CAssignment, n, fill = 0)
# A tibble: 4 × 5
  Reason                             C1    C2    C3    C4
  <chr>                           <dbl> <dbl> <dbl> <dbl>
1 Patient Withdrawal from Domain      2     4     1     0
2 Protocol Deviation                  9     6     8     3
3 Protocol Deviation & Withdrawal     0     0     0     1
4 <NA>                              599   603   274    46
Code
acs_itt_dat %>%
  left_join(pp, by = "StudyPatientID") %>%
  dplyr::count(country, Reason) %>%
  spread(country, n, fill = 0)
# A tibble: 4 × 5
  Reason                             IN    AU    NP    NZ
  <chr>                           <dbl> <dbl> <dbl> <dbl>
1 Patient Withdrawal from Domain      4     3     0     0
2 Protocol Deviation                 19     7     0     0
3 Protocol Deviation & Withdrawal     0     1     0     0
4 <NA>                             1265   115   118    24
Code
save_tex_table(
  make_po_table(
    acs_pp_dat,
    "latex"), 
  "outcomes/primary/acs-pp-outcome-summary")
make_po_table(acs_pp_dat)
n (\%) Low
dose
Intermediate
dose
Low dose
with aspirin
Therapeutic
dose
Overall
Randomised 599 603 274 46 1522
Outcome missing 13 (2.2) 11 (1.8) 3 (1.1) 0 (0.0) 27 (1.8)
Outcome observed 586 (97.8) 592 (98.2) 271 (98.9) 46 (100.0) 1495 (98.2)
Met primary outcome 35 (6.0) 25 (4.2) 19 (7.0) 7 (15.2) 86 (5.8)
Table 18: Primary outcome summary for ACS-PP set.
Data and sampling primary model
mdat <- make_primary_model_data(
  acs_pp_nona_dat,
  c("inelgc3", "agegte60", "ctry"),
  c(10, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- primary_mod$sample(
    data = mdat,
    seed = 81241,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 8,
    parallel_chains = 8,
    adapt_delta = 0.99
))
mdrws <- as_draws_rvars(mfit$draws(c(
  "beta", 
  "gamma_site", "tau_site", 
  "gamma_epoch", "tau_epoch", "y_ppc")))
names(mdrws$beta) <- colnames(mdat$X)

# Transformed samples
epoch_map <- acs_pp_nona_dat %>% dplyr::count(epoch, epoch_lab)
site_map <- acs_pp_nona_dat %>% dplyr::count(site_num, site)
site_facet <- factor(mdat$region_by_site, labels = c("India", "Australia\nNew Zealand", "Nepal"))

Ca <- mdat$ctrA
Cc <- mdat$ctrC

mdrws$Acon <- Ca %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
mdrws$Ccon <- Cc %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
mdrws$Atrt <- mdrws$Acon[-1] - mdrws$Acon[1]
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]

mdrws$OR <- setNames(exp(mdrws$Ctrt), c("Intermediate", "Low with aspirin", "Therapeutic"))
mdrws$OR <- c(
  mdrws$OR, 
  "Ineligible aspirin" = exp(mdrws$beta[7]), 
  "Age \u2265 60" = exp(mdrws$beta[8]),
   setNames(exp(mdrws$beta[9:10]), c("Australia/New Zealand", "Nepal")))
mdrws$ORepoch <- setNames(exp(mdrws$gamma_epoch), epoch_map$epoch_lab)
mdrws$ORsites <- setNames(exp(mdrws$gamma_site), site_map$site)
Code
save_tex_table(
  odds_ratio_summary_table(mdrws$OR, "latex"),
  "outcomes/primary/primary-model-acs-pp-summary-table")
odds_ratio_summary_table(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Intermediate 0.74 (0.42, 1.28) 0.77 (0.22) 0.86
Low with aspirin 0.85 (0.44, 1.56) 0.89 (0.29) 0.70
Therapeutic 2.52 (0.86, 7.13) 2.91 (1.67) 0.05
Ineligible aspirin 2.46 (0.66, 8.05) 2.94 (1.98) 0.09
Age ≥ 60 1.87 (1.14, 3.05) 1.93 (0.49) 0.01
Australia/New Zealand 1.11 (0.24, 4.49) 1.43 (1.18) 0.44
Nepal 1.99 (0.49, 7.13) 2.47 (1.87) 0.16
Table 19: Posterior summaries for model parameters.
Code
p <- plot_epoch_site_terms(mdrws$ORepoch, mdrws$ORsites, site_facet)
pth <- file.path("outputs", "figures", "outcomes", "primary", 
                 "primary-model-acs-pp-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p